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Abstract. The one-point function (i.e., the isotropic flux distribution) is a complementary 
method to (anisotropic) two-point correlations in searches for a gamma-ray dark matter 
annihilation signature. Using analytical models of structure formation and dark matter halo 
properties, we compute the gamma-ray flux distribution due to annihilations in extragalactic 
dark matter halos, as it would be observed by the Fermi Large Area Telescope. Combining 
the central limit theorem and Monte Carlo sampling, we show that the flux distribution takes 
the form of a narrow Gaussian of ‘diffuse’ light, with an ‘unresolved point source’ power-law 
tail as a result of bright halos. We argue that this background due to dark matter constitutes 
an irreducible and significant background component for point-source annihilation searches 
with galaxy clusters and dwarf spheroidal galaxies, modifying the predicted signal-to-noise 
ratio. A study of astrophysical backgrounds to this signal reveals that the shape of the total 
gamma-ray flux distribution is very sensitive to the contribution of a dark matter component, 
allowing us to forecast promising one-point upper limits on the annihilation cross section. We 
show that by using the flux distribution at only one energy bin, one can probe the canonical 
cross section required for explaining the relic density, for dark matter of masses around tens 
of GeV. 
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1 Introduction 


The Large Area Telescope (LAT) onboard the Fermi satellite [1] measured the energy spec¬ 
trum [2] and angular anisotropies [3] of the diffuse extragalactic background of gamma rays. 
Components contributing to this background include blazars [4, 5], star-forming and star- 
burst galaxies [6], and misaligned active galaxies [7]. The combination of these sources gives 
reasonably good fit to the spectral data [5, 8], while the anisotropies are consistent with the 
blazar component alone [3, 9]. Independently of this, dark matter has emerged as the pre¬ 
ferred explanation of many astrophysical and cosmological features through gravity (galactic 
rotation curves, 3> ST&, lensing by galaxy clusters, etc.). If particle dark matter produces 
gamma rays (e.g., by self-annihilation) as in the case of weakly interacting massive particles 
(WIMPs) motivated by popular particle-physics models [10, 11], then it could also contribute 
to this diffuse signal (in some unknown proportion) [12]. Given that the known astrophysical 
sources yield reasonable fit to the spectrum of the gamma-ray background, the dark matter 
component started to be tightly constrained only through the spectral data (e.g., [5, 8, 13]). 

Recently, new analysis techniques beyond the energy spectrum and angular two-point 
correlations were proposed and investigated extensively. Among them is to take cross correla¬ 
tions of gamma-ray data with local galaxy catalogs [14, 15] and matter distribution through 
lensing data [16, 17]. Although recent measurements of the cross correlations [18-21] are 
consistent with the hypothesis of no dark matter signal, they yield tight constraints thereof 
(e.g., [22]). 

Complementary to the studies on these two-point functions, the one-point function 
(i.e. the photon-count or flux distribution) would leverage the isotropic component of the 
diffuse signal. For example, the flux distribution of Milky-Way subhalos has been used to 
constrain particle dark matter properties in light of Fermi unidentified sources [23]. The 
one-point function of Fermi-LAT data has been experimentally fit to a combination of the 
diffuse background and blazar-like sources by Ref. [24], Theoretically, Ref. [25] studied 
the one-point probability density function (PDF) of the gamnra-ray flux due to Galactic 
subhalos, and showed that it features power-law tail at high-flux end (see also Refs. [26, 27]). 
Understanding the one-point PDF for all the relevant sources will be important also for 
possible detections of dwarfs or galaxy clusters with gamma rays (e.g., [28]). 

In this paper, we extend the theoretical framework of Ref. [25] to include the contribution 
of dark matter annihilation in the extragalactic halos. We model the gamma-ray flux from 
the population of dark matter halos using the mass and luminosity functions predicted for the 
structure formation scenario in the Universe with cold dark matter and cosmological constant 
(ACDM). Since there are many such halos, the total flux observed by Fermi is predicted using 
‘large A p statistical tools. By combining the central limit theorem (CLT) in the low-flux 
regime and a Monte Carlo method in the high-flux regime, we find that the differential flux 
distribution has a roughly Gaussian peak, as result of the diffuse emission of a large number 
of very faint sources, but with a power-law high-flux tail due to the rare occurrence of an 
exceptionally bright halo. The all-sky flux in our fiducial model lies around the sensitivity 
limit of the Fermi, well-below the measured gamma-ray background [2], and at roughly the 
same level as the dark matter flux expected from the Fornax galaxy cluster. We find that 
the detectability of a dark matter signature from galaxy clusters over the extragalactic dark 
matter background decreases, when the luminosity boost due to halo substructure increases. 
We also illustrate how to disentangle a dark matter signal from astrophysical backgrounds in 
the presence of the photon shot noise, and forecast the upper limit on the annihilation cross 
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section one might expect to obtain using the one-point function alone (although actually 
performing this analysis with its due rigour is beyond the scope of this theoretical paper). 
Given a fiducial model for the dark matter halo substructure boost, the particle dark matter 
mass, etc., we find a 5cr upper bound roughly a factor of two above the thermal cross section. 

This paper is organised as follows. In Sec. 2, we construct the model of the flux dis¬ 
tribution observed at the Fermi-LAT. In Sec. 3, we detail our specific model choices, and 
in Sec. 4, we present our main results, including a sensitivity analysis of our distribution to 
the model choices, and a probabilistic method for summing the fluxes from a (quite literally) 
astronomically large number 0(1O 22 ) of halos. In Sec. 5, we discuss consequences of this 
study for indirect DM searches, by comparing the predicted distribution to the gamma-ray 
fluxes of galaxy clusters, dwarf spheroidals, and blazars. We conclude the paper in Sec. 6. 


2 Flux probability density function: General formalism 


The goal of the present work is to theoretically predict the PDF P(F), which gives the 
probability of observing a total gamma-ray flux F arising from dark matter annihilation in 
extragalactic halos in a Fermi pixel of a particular size. In this section, we present a formalism 
for constructing P(F) given models for the cosmology, halo properties, and annihilation 
process; we proceed in a completely general manner, postponing specific choices for these 
models until Sec. 3. 

We construct a Bayesian hierarchical model, to predict unknown gamma-ray observables 
from well-constrained ACDM parameters and fitted models of N-body simulations. In the 
hierarchical Bayesian approach, uncertainties in the parameters of probability distributions 
are modelled with their own distributions (and recursively). This allows us to systematically 
combine the uncertainties on physics at widely differing scales, and thereby to perform a 
sensitivity analysis of our model (Sec. 4.1). 


2.1 PDF for the flux from individual halos 

Throughout the paper, F represents the differential flux, i.e., a number of photons received 
per unit area, unit time, and unit energy range [F(E) = d 3 Nj/dAdtdE\. The PDF P(F ) 
for observing a total differential flux F from all of the halos in a pixel depends on the PDF 
P\(F) for observing F from any individual halo. 1 We thus proceed by first deriving the latter 
quantity. 

Because the differential flux F from an individual halo is completely determined by its 
rest-frame differential luminosity L = dpN^/dtdE and its redshift z, we can write 


Pi(F) = jdL dz P(F\L, z)P(L, z) 


= J dLdz5[F - F(L,z)]P(L\z)P(z) 

Here, the usual relation for the differential flux, 


( 2 . 1 ) 


( 2 . 2 ) 


1 Throughout this paper, we denote probability distributions by P(- ■ ■) and distinguish them using the 
random variables that they describe, along with subscripts if necessary. Conditional and parameterised dis¬ 
tributions are denoted as P(-|-)- Exceptions to this convention are Poisson and Normal distributions, denoted 
V (-I-) and Q |, •) respectively. 
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depends on the luminosity distance d,£,(z ) and the pair-production optical depth t(E,z) 
for gamma-ray photons, and also accounts for the redshift of photons emitted with rest- 
frame energies E( 1 + z) to observed energies E. We can interpret P(L\z) = dN/dL(z ) as 
the redshift-dependent halo differential-luminosity function. 2 Assuming that the halos are 
isotropically distributed across the Universe, the number of halos at redshift z is proportional 
to the comoving volume 5V(z ) of the corresponding redshift slice Sz, therefore we also have 
P(z) = dN/dz oc dV/dz. 

Alternatively, we can rewrite Eq. (2.1) in terms of the halo mass M to obtain 

Pi(F) = JdLdMdz6[F - F(L,z)\P(L\M,z)P(M\z)P(z), (2.3) 

where we can similarly interpret P{M\z) = dN/dM(z) as the redshift-dependent halo mass 
function. 

In principle, the distribution P(L\M, z) in Eq. (2.3) captures the scatter in the relation 
between the differential luminosity and the mass of a halo, which also depends on redshift. 
This is because the halo luminosity is determined not only by the properties of the dark 
matter particle and the details of the annihilation process, but also by the density profile p of 
the halo, which usually shows scatter for any given M. The halo profiles can be completely 
characterised by some parameters Oh (such as p s ,r s , r v i r ... in the case of the NFW profile [29]) 
so that (for any given particle dark matter model) we have L = L(Oh)- If we further assume 
that the distribution of halo profiles can be described by a halo model that gives P(Oh\M, z), 
we can write 

P(L\M,z) = J dO h P(L\O h )P(O h \M,z) 

= Jd0 h 6[L-L(0 h )]P(O h \M,z). (2.4) 

We can then use this expression to simplify Eq. (2.1), giving 

Pi(F) = JdMdzd0 h 5[F-F(0 h ,z)}P(0 h \M,z)P{M\z)P(z), (2.5) 


where the flux relation is now written in terms of Oh- 

In order to make the numerical calculation of Eq. (2.5) more tractable, we shall neglect 
the scatter in the distribution P(0h\M, z) in this work. That is, we take the distribution of 
the halo-profile parameters 0f L = {Oh, i, ■ ■ •, 0h, n } to be given by 

n 

P(O h \M, z) = H S[0 h ,i - O h ,i{M, z )}, (2.6) 

i =1 


where the functions Oh,i(M, z) give the mean values for the parameters. With this assumption, 
we can perform the integrals over Oh and M in Eq. (2.5), leaving only an integral over z: 


Pi(F) 



dF _1 dN dV 
~8M dM~dz' 


(2.7) 


"Here and elsewhere, equalities between probability densities and number densities is meant modulo a 
normalisation. 
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Here, functions of M in the integrand are evaluated at the value of M defined implicitly by 
the flux relation F = F(9h(M, z ), z) for the corresponding values of F and z. 

Eq. (2.7) then gives the PDF Pi(F) for the differential flux from individual halos. To 
reiterate, this is given in terms of (i) the cosmological model, which affects the differential 
flux F(Oh,z), the volume P(z), and the mass function P(M\z) w , (ii) the halo model, which 
gives the mean values z) of the halo-profile parameters; and (iii) the optical depth and 

the details of the annihilation process, which affect the normalisation of F(9h,z). In Sec. 3, 
informed by observations and simulations, we shall make specific, fiducial choices for these 
quantities and calculate the resulting P\ (F). However, before doing so, we shall complete 
our discussion of our general formalism by reviewing how P\ (F) can be used to find P(F), 
the PDF of the total flux F from all of the halos in a pixel. 

2.2 PDF for the total flux 

We assume that all dark matter gamma-ray sources may be treated as point sources (see also 
Sec. 3.1.1 and Sec. 5.4 below), allowing us to equate the differential flux F and the differential 
flux per pixel. The dark matter differential flux F arriving at any given pixel of the Fermi 
sky map, is the summed flux F = JT Fi of any number of individual halo point sources [30], 
where each differential flux F t is an independent and identically distributed (i.i.d.) random 
variable with the distribution Pi(F). The distribution of a sum of random variables is the 
convolution of all the original distributions [31]; Since the F t are i.i.d., the distribution of the 
total differential flux per pixel is the autoconvolution [25] 

P k (F) = P 1 (F) *Pi(F) * ■ ■ ■ * Pi{F) = (Pi)* fc , (2.8) 

where k is the number of halos contributing to this flux. Since furthermore we do not know 
how many halos are thus stacked in a pixel, the number k of fluxes in the sum is itself a 
random variable. If we assume this number k of halos per pixel is Poisson-distributed over 
the sky with some mean N', we can model the total differential flux per pixel as 

P(F) = f dN'P(N') J2 P ( k\N ') P k (F) , (2.9) 

•' k 

where the uncertainties in k and N' are marginalised away. Since the numbers k and N' of 
extragalactic halos are very large, both P(N') and V ( k\N') are thin enough to be approx¬ 
imated by delta functions, so that P(F) = Pk^N'{P)- Thus, the only additional physical 
input required to compute P(F) from Pi(F) is N 1 , which is discussed below. 


3 Model inputs of the flux probability density function 


3.1 Cosmological inputs 

3.1.1 Number of halos per Fermi pixel 

Cosmology directly determines the redshift distribution of halos (via isotropy P(z) = dN/dz ) 
and their mass distribution (via the gravitational collapse of inhomogeneities that yields 
dN/dM ). The normalisation of these number densities clearly corresponds to the total num¬ 
ber N of halos in the Universe. The number of halos per Fermi pixel is then 


N' = 


= 

4-7T 


fi 


pix 


47T 


dN dV 
dM dz 


dMdz , 


(3.1) 
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where O p i x is the pixel size (expressed in units of a solid angle) and the mass function dN/dM 
is described in Sec. 3.1.2. The parameters of the ACDM cosmology relevant to our model 
are (Planck+WMAP from [32]) 

{n b h 2 , n c h 2 , n A , H 0 , a 8 ,n s } = { 0 . 0221 , 0 . 120 , 0.685, 67.3,0.83,0.96} . ( 3 . 2 ) 


Furthermore, the integration limits for Eq. (3.1) (and also Eq. (2.7)) are chosen as follows: 
we assume that our dark matter candidate forms structures down to lO _ 6 A /0 [33-37], and 
allow for virialised dark matter structures up to 10 17 Mq, 100 times more massive than galaxy 
clusters. We assume that we can measure luminosity from structures that form between z = 5 
and 10~ 5 , the latter of which corresponds to a distance of roughly 45 kpc, well outside of the 
baryonic content of the Milky Way in any direction. 

We model the energy-dependent angular resolution of Fermi-LAT [1] as follows: 


0(E) 


0.8° (E/ GeV ) -0 ' 68 for 0.04 < E/GeV < 20 
0.1° for 20 < E /GeV 


(3.3) 


such that 9 ss 0.8° at our fiducial observing energy (justified in Sec. 4.1) of E = 1 GeV. 

This is slightly larger than 0.6° quoted in [1] , which is valid for normally incident photons 
only. Since we adopt this angular resolution as a size of each pixel, fl p i x ~ it9 2 (E), we entirely 
neglect the instrumental point-spread function. This is also justified because the number of 
sources per pixel is found very large ( N' ~ 7 x 10 21 ) and the flux is diffuse. 

The largest possible mass of a point source at a given redshift M-^xt(z) may be deter¬ 
mined from the critical virial radius i?Ext ~ Al(z) tan(0/2) that fits in a pixel. We can use 
this as an integration limit in Eq. (3.1) to cross-check our assumption (Sec. 2.2) that all dark 
matter gamma-ray sources are point-like. We find less than 0.28 extended sources per pixel. 
This is much more than twenty orders of magnitude smaller than the total number of sources, 
but still represents a non-negligible absolute amount given the large number of pixels. We 
discuss this issue more thoroughly in 5.4. 


3.1.2 Halo mass function 


The halo mass function, first addressed heuristically by Press and Schechter [38] and subse¬ 
quently formalised in, e.g., Ref. [39], is computed as 


dn 

dM 


Mm— 

M jk 'dM 



(3.4) 


where 5 C = 1.69 /D(z) is the (linear) critical overdensity, D(z) is the linear growth factor, and 
o(M) is the mis deviation of primordial density fluctuations, smoothed to scale M [39]. The 
functional form of a(M) (required to calculate du/dM ) is determined from the literature 
[40] with normalisation set by the cosmological parameter a 8 [32], The function f(v) is 
derived from the excursions of these density fluctuations above a ‘barrier’ [39] that encodes 
the physics of halo collapse (including <5 C ). For an approachable presentation of the formalism, 
see Ref. [41]. 

In addition to ellipsoidal collapse, our fiducial mass function incorporates a virialised 
halo’s angular momentum and the cosmological constant into its barrier 8 C . It has a self¬ 
similar f(u) well-fit by the following function (Eq. (163) in [42]): 


vf(v) oc 



— 0A019au 


0.5526 

1 + 7 -+ 


0.02 


1 2 > 


(az/) 


0.585 


(H 


0.4 
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where a = 0.707. The resulting mass function is similar to the more common Sheth-Tormen 
parameterization [43, 44] (within the resolution of existing simulations), demonstrating that 
the extra physics of our barrier have very little effect on the high-mass end of the mass 
function. However, we expect the cosmological constant to delay the formation of large-scale 
structure, leaving us with a larger proportion of halos at high redshift that are smaller than 
current simulations can resolve. We find this gives roughly three times more flux than if we 
had used a mass function for which the cosmological constant is ignored. 

3.2 Halo model 

The differential luminosity from annihilation in a dark matter halo of mass M is given by the 
product of a particle physics term and an astrophysical J-factor, i.e. the line-of-sight integral 
of the dark matter density squared, boosted by the annihilations in halo substructures. The 
dark matter density can be parameterised using the same profiles that fit N-body simulations 
well. For an NFW profile, the J factor of a point-like source can be analytically recast as 

J = (1 + B) a(c vir )p s M vir , (3.5) 

in terms of the substructure boost B , the virial concentration c v ; r , the scale density p s of the 
profile, and the analytic function (e.g., [45]) 

“ (c) = (‘ _ (TT^f) ( ln(1 +c) ~ rb) ■ 

The concentration parameter c of an NFW halo is related to the background density at 
the time that the halo forms: small mass halos are more concentrated than high-mass halos 
because they form earlier (hierarchical halo formation). The concentration is also the link 
between scale parameters ( p s ,r s ) of the halo profile and the halo’s mass content [46]. 3 * * 

The presence of halo substructures enhances the luminosity of the halo as a whole. 
Furthermore, halo substructure is expected to be denser than a host halo of the same mass 
(e.g., [45]). This higher density entails a larger number of annihilations, further enhancing the 
luminosity. The boost factor B parameterises this substructure luminosity as a proportion 
of the host luminosity. Since the J factor increases as the density squared, substructure is 
expected to contribute between twice to twenty times as much luminosity as host structures 
of mass 10 6 < M/M 0 < 10 13 [48]. 

There are a number of N-body fitted models for NFW concentration and substructure 
boost to choose from in the literature. The models for B used in this study are listed in 
Table 1. 

The optimistic model is based on fit to the numerical simulations [49], which is well 
motivated for mass scales larger than the current resolution limit. For smaller scales, on the 
other hand, it heavily relies on validity of the phenomenological extrapolation. Reference [48] 
pointed out that such a power-law extrapolation was unphysical, and came up with more 
physically motivated model by adoping the flattening of the concentration-mass relation 
found for held halos. 

3 We convert between the virial and (• • ■ Roo conventions found in the literature using the prescription of 

Ref. [ 47 ]. This allows us to interchangeably convert between c v n and C200 at a given mass, and to convert 

between M v i r and M200 given a choice of concentration-mass relation. 
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Boost models 

[ref] 

Formula 

No boost 


B = 0 

Fiducial 

[48] 

Eq. (3.6) 

Optimistic 

[49] 

B = 1.6 X 1O- 3 (M 2O o/M 0 ) 0 - 39 


Table 1: Substructure boost models considered in this study, ordered from least to most 
optimistic in terms of annihilation signal detection prospects. The nonzero models are derived 
from fits to N-body simulations. 


The boost factor for our fiducial model is given as [48]: 

5 

logioB = £^°gio(M 2 ° o r, (3.6) 

i=0 

B, = {-0.442,0.0796, -0.0025,4.77 x 10" 6 , 4.77 x 10~ 6 , -9.69 x 1(T 8 } . (3.7) 

The underlying concentration model is very close to other N-body-motivated models [50, 51]. 
The third model represents the most conservative (but unexpected) situation in which there 
is no substructure boost. In this study, we consider all the boost models in order to bracket 
uncertainties on subhalos, although the fiducial model is preferred over the others. 


3.3 The gamma-ray model 

3.3.1 Dark matter annihilation model 

The particle physics component K(E) = (av)(dN/dE )/of this model ( L(E ) = JK(E)) 
is assumed here to be a standard WIMP model, with annihilations into gamma rays via bb. 
(av) is taken to be the thermal cross section 3 x 10~ 26 cm 3 s _1 , and the WIMP mass is taken 
to be m x = 85 GeV. Our parameterisation of the photon number per energy per annihilating 
particle is [52]: 

dN 0.42 exp(—8x) E 

dE m x (x 1 - 5 + 0.00014) ’ X m x ' 

With the values quoted above, we have K = 1.3 x 10 -38 cm 3 s _1 MeV - ' 3 at 1 GeV. Since 
our method is independent of any specific particle physics model, and depends only linearly 
on K, it is trivial to rescale this study’s results to accommodate other particle physics models 
(as demonstrated in Sec. 5.3). 


3.3.2 Gamma-ray optical depth 

By restricting our analysis to small enough redshifts and energies, we do not need to consider 
photoionisation or pair production [53]. We can then use a very rough parameterisation [52] 
of gamma-ray absorbtion: 


g —t(E, 2) = exp 


z f E 
T3 V10GeV 


0 . 8 ' 


(3.9) 


That this is very simplistic does not matter too much, since it is not a very important effect 
at low energies and redshifts anyway — as can be justified quantitatively by a sensitivity 
analysis. 
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4 Results 


In this section, we calculate P\ (F) and P{F) given the physical inputs of Sec. 3. We find that 
P\{F ) is roughly power-law like with a log-slope of approximately —2 over a broad dynamic 
range, with some model-dependent features to which P{F) is sensitive. More importantly, 
we find that P(F ) takes the shape of a Gaussian at low flux, connecting smoothly with a 
power-law tail at high flux. This power-law tail, which has the same slope as Pi(T’), is the 
regime where the flux is dominated by a singularly bright source. 

4.1 Sensitivity analysis of Pi(F) 

The probability [Eq. (2.5)] that any single halo produces a differential flux observed with a 
value F is then given by marginalising away the uncertainty in its mass and redshift. This 
marginalisation is clearly sensitive to the model choices we have presented above. We now 
discuss the sensitivity of P\ (F) to the model choices. 

4.1.1 Dominant effects 

We find that the two most significant effects on P\ (F) (from amongst the effects we consid¬ 
ered; see below) are the choice of observing energy at the Fermi-LAT, and the substructure 
boost model B(M). These are represented in Fig. 1, which illustrates not only that Pi(F) 
is more complex than just a power law P\ (F) oc F~ 2 (even in the absence of a substructure 
boost), but also that the substructure drastically influences the shape of P\(F). 

The choice of the energy at which we study our differential flux represents a tradeoff 
between the amount of flux we expect to observe (which decreases with energy; Fig. 1), and 
the instrumental capacity to actually observe it (which increases with energy; Ref. [1]). The 
observing energy of 1 GeV adopted in this study represents the lowest energy at which we 
can leverage Fermi’s best angular resolution and effective area. The instrument’s energy 
resolution of A E/E = 9% at 1 GeV is optimum for normally incident photons [1]. Any 
systematic error due to misidentified energies is therefore minimised by this choice of energy 
(although we do not account for this instrumental effect in our analysis). 

Consequently, the effect of a gamma-ray optical depth (which predominantly affects 
photons from distant sources) only affects the low-flux region of P\(F). This simply reflects 
the fact that at equal luminosity, distant sources produce less flux. At our low energy of 
1 GeV, the attenuation factor of a source at z = 3 is e~ T ~ 0.86: the net effect is definitely 
smaller than a few percent after marginalising over the nearby halos, for which the attenuation 
is truly negligible. A better parameterisation than [52] was therefore not deemed necessary 
for this exploratory analysis. 

Since our differential flux is roughly F ~ L/d? Ll there is an extent to which ‘bright’ 
halos are massive and nearby, while ‘faint’ halos are light and distant. Therefore the largest 
modeling factor in the high-flux tail is the choice of boost model. This is reflected in the right 
panel of Fig. 1: the optimistic model [49] gives very large high-flux tail, in comparison to the 
fiducial model. The mass scale at which the fiducial and optimistic boost models intersect 
(~ 10 9 Mq [48]) is present in the flux distribution. We find that Pi(F) is less sensitive to 
other modeling choices. 

4.1.2 Subleading effects 

We have used various analytical models for the power spectrum transfer functions T(k ) 
[54, 55], and for c v - u fits to N-body simulations [48, 50, 51]. The resulting distributions P\{F) 
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Figure 1: The flux/intensity PDF for a single dark matter source, P\(F), with its depen¬ 
dence on photon energy (left) and boost models (right). In the right panel, the blue, black, 
and red curves represent respectively the pessimistic, fiducial, and optimistic models of the 
subhalo boost. The choice of boost model clearly and significantly affects the functional form 
of the one-point function. The log-slopes of the fiducial model are offset (black dashed) and 
quantified for convenience. The flux F and intensity I of the gamma-ray background are 
related via the instrument’s pixel size: F = /fl p i x , where D p i x ~ 5.8 x 1CD 4 sr for E = 1 GeV 
photons. 


were found to be robust to changes in these inputs. Indeed, the function a(c) in the J 
factor [Eq. (3.5)] is very smoothly decreasing for c > 1, so the choice of concentration model 
does not influence the final result much. Similarly, a(M) only changes by about an order 
of magnitude over many orders of magnitude of halo masses, so varying models of T{k) (or 
even using an N-body fit [56]) does not significantly change P\(F) either. 

Using such a simple, ‘self-similar’ mass function, is not without shortcomings: these 
Markov-process models do not follow individual halo histories, and cannot account for dy¬ 
namical effects such as dynamical friction or tidal stripping. This, amongst other concerns, 
underlines the danger of our uncontrolled extrapolation of the mass function down to halo 
masses 10~ 6 Mq. Nevertheless, mass function fits to N-body simulations are often very good; 
the fits even favour ellipsoidal collapse models [43, 44] (which we adopt in this study) over 
the spherical collapse models of the seminal papers. The ellipsoidal collapse model is accu¬ 
rate to a few % over the large halo mass range (while spherical collapse [38] underestimates 
the number of FOF halos [57]) and gives an excellent fit in the range of 10 5 -10 9 Mq [58]. 
Ellipsoidal collapse is also suitable for halos as small as 10 3 Mq and as early as z = 15 [59]. 
We studied dependence of the mass function on A, in addition to ellipsoidal collapse [42], 
Without significantly changing P\ (F).we find that the delayed structure formation gives a 
roughly three times larger total flux (by increasing the number of halos N) than if we had 
used a mass function for which the cosmological constant is ignored [43, 44], 

We will need, when computing the total P(F), to compute the first few moments of 
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P\(F) at an intermediate step: 


P] (F) = J dE , Pi(E , )E', 

(4.1) 

Pl ( F ) = J dF'P 1 {F')(F'- E Pl(F) ) 2 . 

(4.2) 


After multiplying Eq. (4.1) by the mean number of halos [Eq. (3.1)] and dividing by the 
pixel size, one obtains the mean intensity of the gamma-ray background from dark matter 
annihilation [60]. Similarly the variance [Eq. (4.2)] is related to the angular power spectrum 
after similar corrections [61]. 

Since Eq. (2.7) entails an integration over redshift for each F, and since we find Pi(F) 
to be relatively smooth, we calculate 250 logarithmically equidistant points over the ~40 
orders of magnitude supporting the distribution. In order to obtain robust estimate of the 
moments, we further sample 250 points within the four orders of magnitude nearest to the 
maximum estimated from the low-resolution sampling. 

Sources of uncertainty not studied include the cosmological parameters (3.2), instru¬ 
mental effects related to Fermi-LAT, and the assumed NFW profile with exactly determined 
parameters 0(M,z). Of these, the largest quantifiable uncertainty is the scatter about the 
concentration C 20 O) estimated at 15% [51]. 

4.2 Computing E(E) 

4.2.1 Monte Carlo method combined with the central limit theorem 

Once P\ (F) and the number of halos per pixel N' is specified, the calculation of E(E) requires 
no additional physical assumptions. However, the large number of halos per pixel k ~ N' 
makes any exact calculation of the autoconvolution [Eq. (2.8)] prohibitively expensive, even 
using Fourier methods [25, 30, 62], One might attempt to use the central limit theorem 
(CLT) to approximate Pk(F) by a Gaussian. But, although the CLT guarantees convergence 
in distribution for k —> 00 , at finite k there will be deviations from a Normal distribution, 
particularly in the tails [63]. This is especially true for power-law-like distributions such as 
Pi(F), since the stable distributions of sums of power laws may be non-Gaussian [31]. See 
also Appendix A. 

For the purpose of solving the autoconvolution [Eq. (2.8)], we can choose to split P\ (E) 
into two physically interesting contributions: the multitude of low-mass, faint halos con¬ 
tributing less than some cutoff flux E*; and the rare, large point sources brighter than E*. 
This scheme is illustrated in Fig. 2. The former will contribute an isotropic background 
with Gaussian statistics, while the latter will add a power-law contribution to the high-flux 
tail, which smoothly matches with P\{F) at high fluxes (where flux is dominated by a sin¬ 
gle source). The noise at high flux in our Monte Carlo generated E(E) can therefore be 
confidently ignored as a numerical artefact (cf. Fig. 10 in Appendix B). 

Despite the large value of k , our cutoff E* may be chosen such that only a few thousand 
of these high-flux sources remain in each Fermi pixel. We can model the contribution of these 
sources by Monte Carlo, drawing number of these rare sources from a Poisson distribution 
and their flux from Ei(E). We obtain finally the flux from k halos as the sum of these two 
contributions; The flux distribution from k sources is given by the following convolution: 

Pk{F) = (E) * P*>F. (E) . (4.3) 
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Figure 2: Schematic of the F ,* cutoff of P\ (F) into high/low flux. In our computation of 
the full P(F), the central limit theorem is used to combine the fluxes from the many sources 
fainter than _F*, that follow a distribution P^ <F * (F). Monte Carlo is used above this cutoff 
to combine the halo fluxes drawn from Pf >F *(F). 

We note that when taking the convolution of our Monte Carlo result with our faint- 
halo Gaussian, we care mostly about the peak of the Gaussian (since the peak is the largest 
contribution to the convolution integral). The CLT thus offers a suitable approximation of the 
Pf <F *(F) autoconvolution for this purpose, despite deviations from a Normal distribution 
in tails. In fact, for practical purposes only the peak matters and the convolution, Eq. (4.3), 
is a trivial shift of the Monte Carlo result to higher fluxes. 

For a more detailed derivation of these results, see Appendix B. 

4.2.2 Flux distribution and instrumental sensitivity 

In Fig. 3, we show the flux PDF P(F) for the three subhalo boost models. As discussed above, 
the distributions are well represented by the ‘diffuse’ component of nearly Gaussian with a 
power-law tail at high-flux regime. The mean of this distribution [Eq. (4.1)] corresponds 
to ~10~ 12 cm~ 2 s -1 sr -1 MeV -1 of the diffuse gamma-ray background, with an order-of- 
magnitude level dependence on the model (see Table 2). This is the value used as the mean 
intensity in the literature, in order to constrain the dark matter annihilation cross section 
from the comparison with the spectral data [64]. However, our PDF analysis shows that the 
distribution is skewed, such that the mean is not the most likely value to be observed in any 
given pixel: the mode is typically lower than the mean, by a boost-dependent factor of the 
order of a few percent, again summarised in Table 2. If one instead uses these most likely 
values, then the upper limits on annihilation cross-section will accordingly remain relatively 
stable: existing upper limits are thereby relatively immune to the non-Gaussianity of the 
P(F) tails. However, for accurate results, one has to perform the data analysis by taking 
into account the full shape of P(F). 

Before contrasting our dark matter signature to known and well-observed astrophysical 
sources such as galaxy clusters and blazars in the next section, we briefly touch on whether 
Fermi is sensitive enough to see it at all. A flux of a single, GeV photons per pixel, over 
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Figure 3: The flux PDF P(F ) per pixel. The blue, black, and red curves represent respec¬ 
tively the pessimistic, fiducial, and optimistic models of the subhalo boost. Instrumental 
responses of Fermi-LAT on detecting P(F) are schematically shown. Vertical lines represent 
a flux corresponding to a single, one GeV photon per pixel, over the course of a mission of 
duration 5 (10) years. The Horizontal line schematises the angular resolution limit [Eq. (4.4)] 
at 1 GeV. 


Boost model 

Mean 

Most Likely 

Difference 

Ratio 

EGB fraction 

No boost 

1.0 

1.0 

0.0 

x 1.0 

0.2% 

Fiducial 

3.68 

3.52 

0.16 

x 1.05 

0.6% 

Optimistic 

15.2 

11.9 

3.3 

xl.3 

2.5% 


Table 2: Mean and most likely extragalactic dark matter annihilation intensities as a func¬ 
tion of the substructure boost model, in units of 10~ 12 cnr 2 s -1 sr _1 MeV -1 . The difference 
between these quantities is percent-level, securing existing constraints on particle dark mat¬ 
ter properties against the non-Gaussianity of P{F). We also provide, for interest, the value 
of this mean contribution of the DM as a fraction of the unresolved EGB at 1 GeV [65]. 

a 5 year mission with LAT’s effective area of 0.9 m 2 and a field of view of 1/5 of the sky, 
corresponds to a differential flux of 6 x 10~ 12 cm~ 2 s -1 sr _1 MeV -1 . The bulk of the one- 
point function P(F), with its peak of 3.5 X 10~ 12 cnr 2 s -1 sr -1 MeV -1 , lies just below this 
sensitivity limit. See the vertical dashed lines in Fig. 3, for the sensitivity curves for 5-year 
and 10-year Fermi exposure. 

There will be a small fraction of the pixels that register photons from the high-flux tail. 
Since the high-flux power-law tail is characterised by Pi(F) (see Fig. 10 in Appendix B), 
this is to some extent equivalent to computing “the probability of seeing a dark matter point 
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source.” The finite angular resolution of the LAT limits the number of pixels in which we 
can look for these bright outliers: we argue in Appendix C that a maximum flux F max , the 
brightest flux one might expect in any pixel, is given by 

F max P(F max ) « , (4.4) 

-^pix 

where -/V p i x = 4it/(it0 2 ) = 2.18 x 10 5 is the number of pixels. This equation simply states 
that we are unlikely to probe the flux regime of P{F ), where odds are worse than 1 /A r p ; x . In 
Fig. 3, we show F max P(F max ) as a horizontal dashed line, and this confirms that one is able 
to probe this high-flux tail with Fermi-LAT’s very good angular resolution. 

Since the fiducial model lies on the edge of detectability with Fermi, it would not 
substantially contribute to the observed extragalactic background with the fiducial choices 
of dark matter parameters; this is consistent with what is found in the literature (e.g., [60]). 
Of course, the particle physics parameter space [Eq. (3.8)] may still allow such contributions 
for e.g. larger values of ( ov } than the thermal cross section: in Sec. 5.3 we will vary (crv) 
slightly above the fiducial value, and therefore expect a flux of a few photons per pixel over 
the course of the 5-year mission. Similarly, with an optimistic boost model, the dark matter 
may contribute a few photons despite a using the canonical particle physics model (Fig. 3). 


5 Discussion 


5.1 Searches for clusters and dwarf spheroidal galaxies 

This diffuse emission due to dark matter annihilation characterised by P(F) acts as a back¬ 
ground to dark matter point source searches in galaxy clusters and dwarf spheroidal galaxies. 
The situation is summarised in Fig. 4: The flux from promising candidate sources is su¬ 
perposed onto the extragalactic dark matter one-point function, for the three boost models 
presented in Table 1. A discussion of astrophysical backgrounds is postponed until Sec. 5.2.1: 
indeed, it is unlikely to find by coincidence a bright blazar in the same pixel as a drawf galaxy. 


5.1.1 Clusters of galaxies 


Our candidate sources in this category are the Fornax and Coma clusters [66, 67]. These 
targets are among those with the most promising J-factors, so that they represent the most 
optimistic prospects of an annihilation signal detection (or alternatively the strongest con¬ 
straints on such a signal). 

The substructure of a halo lies predominantly outside the scale radius r s = r v ir/c v ir [49, 
68]. In the absence of substructure boost, in contrast, almost all the luminosity is concen¬ 
trated within this scale radius and our clusters are well approximated as point sources for 
Fermi-LAT: 


9b= o = arctan 



< e 


pix j 


(5.1) 


where the fiducial concentration model c v i T (M,z) is used, and cIa is the angular diameter 
distance. 

If the boost is significant, then substructure outside the scale radius of these clusters 
contributes to the luminosity: light comes from the entire virial radius, and the flux from 
the halo (treated as an extended source) is diluted between many adjacent pixels. The boost 
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Figure 4: Brightnesses of promising clusters and dwarf spheroidal galaxies superposed on the 
extragalactic dark matter annihilation gamma-ray background. The color code is the same 
as for previous figures. We assume that dwarf spheroidals have no substructure boost. The 
fiducial model does not favour indirect searches with clusters. The inversion of predictions for 
Coma and Fornax between top and bottom panels accounts for source extension, as explained 
in the main text. 


factor gives the fractions of luminosity from the host and its substructure, so a rough estimate 
of the cluster radius we should convert into an angular extension is 


R = 


f s T Br v j r 
1 + B 


(5.2) 


We clearly recover R ~ r s when the contribution from substructure is negligible (B <C 1), 
and R ~ r v - ir when B 1. The angle arctan {R/d^) then determines the number of pixels 
over which the flux is averaged into an intensity. This corresponds to flux dilutions over 
roughly 10 pixels for Coma and 60 pixels for Fornax, explaining why intensities from Coma 
and Fornax appear inverted in lower two panels of Fig. 4: the total flux increases when 
considering substructure, but flux per solid angle decreases more for Fornax than for Coma. 
The fact that the intensity from Fornax appears to decrease from the top panel to the lower 
panels of Fig. 4 is then just a manifestation of the difference between seeing Fornax as a 
point source in the top panel or as an extended source in the lower panels. 

For the optimistic boost model, Coma stands out in the tail of P(F), while Fornax 
is only barely more visible than if it (pessimistically) had no substructure. Although our 
treatment of source extension is somewhat naive, the diffuse gamma-ray background would 
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be a limiting factor in cluster analysis for the fiducial boost model, even if the effects of 
extension were favourably revised by a factor of three or four (see middle panel of Fig. 4). 

The intrinsically poor signal-to-noise in cluster searches is also independent of the anni¬ 
hilation cross section or mass: changing these particle physics parameters would not change 
the signal-to-noise ratio, since both the target cluster and the gamma-ray background are 
rescaled by the same factor. 


5.1.2 Dwarf spheroidal galaxies 

Dwarf spheroidal galaxies are another strategic choice for gamma-ray point-source searches, 
to which dark matter annihilations in the Milky Way substructures constitute a known 
background [28]. We now discuss the background due to extragalactic sources. 

We consider Draco and Segue 1, with J-factors from Ref. [69] since we do not expect our 
virialised halo model (Sec. 3.2) to apply to them. Again, these sources have larger J-factors 
than other known dwarfs, and should thereby set the strongest constraints on non-detection. 
These sources lie well above the isotropic background component, even if we assume no 
relevant substructure boost for the dwarfs (B = 0). Consequently, the more substructure 
boost there is in extragalactic dark matter halos, the worse the signal-to-noise for dwarfs will 
be. Since we know the full distribution P(F ) for the extragalactic background, the p -value of 
an excess signal at 1 GeV (due to a dwarf spheroidal) can readily be estimated. Even though 
the mean intensities yield poor signal-to-noise ratios, a pixel as bright as Segue 1 would be 
relatively uncommon (though not absent from the Fermi skyrnap, see Fig. 3). 

5.2 Astrophysical backgrounds: Blazars and other components 

In this section we consider a number of known astrophysical backgrounds that would mask 
the signal of our fiducial model. 


5.2.1 The unresolved blazar flux distribution 

Following the parameterisation of the blazar source count from Ref. [4], we assume the 
following power-law for P\ (F) of unresolved blazars: 


d 2 N 

dSdQ 


= 1.4 x 1(T 7 



-1.64 

cm 2 s deg -2 , 


Pi{F) = 2.5 x 1(T 10 


F 


cm -2 s -1 MeV” 1 


-1.64 

cm 2 s MeV, 


(5.3) 

(5.4) 


where S is the gamma-ray flux integrated above 100 MeV. We extrapolate this relation down 
to a lower flux limit of S > 0.36 x 10” 10 cnr 2 s _1 , and perturb around this fiducial value 1 
to test the model’s sensitivity to this extrapolation (see Table 3). We also enforce an upper 
limit of S < 2 x 10” 8 cm' 2 s _1 , in order to maintain a low blazar detection efficiency [4] 
without altering the power-law form of the blazar Pi(F). Assuming an JV 2 ' 4 spectrum, the 
corresponding range in differential flux at 1 GeV is then 2.0 x 10” 15 < F/{cvoT 2 s” 1 MeV -1 ) < 
1.1 x 10” 12 . 

Using d 2 N/dSdQ, we And 2.0 faint blazars per pixel, and draw the number of sources 
per pixel from a Poisson distribution in order to predict P(F) by Monte Carlo simulation. In 
practice, we split the blazar P(F ) into two contributions: the delta-function of zero sources 

4 This value is ten times smaller than the ‘medium band’ from Ref. [4], corresponding to the faintest observed 
source out of the mix of FSRQs and BL Lacs that contributes at 1 GeV. 
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Intensity F/£l p ; x [cm 2 s 1 sr 1 MeV x ] 

Figure 5: One-point function P{F ) for the three dark matter models (with boosts color- 
coded as previously), alongside the P(F ) of the diffuse contribution of blazars (green). The 
dashed red band represents the measurement of the unresolved EGB from the Fermi data at 
1 GeV [65], while the dashed green line is the mean of the blazar PDF. 


Smin 

blz/pix 

(F)/Qpi x 

(% EGB) 

0.72 

1.25 

0.89 

(14.8%) 

0.36 

1.97 

0.921 

(15.3%) 

0.18 

3.08 

0.946 

(15.8%) 


Table 3: Sensitivity summary of our unresolved blazar model (with fiducial values in the 
central row). The first two columns pertain to Eqn. (5.3): S m j n (in units of 10~ 10 cm~ 2 s -1 ) 
is the lowest flux to which we extrapolate the source count distribution, from which a number 
of (faint, unresolved) blazars per pixel may be derived. The next two columns summarise 
the corresponding mean intensity of the blazar P(F ), both as an absolute value in units of 
ICE 10 cm~ 2 s' 1 sr -1 MeV -1 , and as a proportion of the unresolved EGB [65]. 

per pixel (treated analytically) and the contribution of more than zero sources per pixel 
(Monte Carlo): 

OO OO 

PBlazar(^) = J2 V W 2 ' 0 ) P ^ F ) = ^^(F) + ^ V (fc|2.0) P k (F) . (5.5) 

k =0 k =1 

The (nonzero) blazar flux distribution is plotted in Fig. 5. Like the flux distribution for the 
dark matter, it has a non-negligible skew and approximates the single-source distribution at 
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high fluxes. i 3 Biazar(^ ? ) also follows the single-source powerlaw between the minimal intensity 
of a single source I m i n = F m i n /Q p i x and the minimal intensity of two sources, with a discon¬ 
tinuous derivative at this threshold. We see from Fig. 5 and Table 3 that unresolved blazars 
contribute roughly 15% of the unresolved EGB at 1 GeV, reproducing Ref. [4]: even a single 
blazar is at least as bright as the entire dark matter contribution of our fiducial model, and 
the mean blazar flux is between one and two orders of magnitude brighter than the dark 
matter, depending on the boost model. Fortunately, using the entire distributions (instead 
of just their means) allows the thin, peaked dark matter to be statistically extracted from 
the broad, powerlawlike blazars (see Sec. 5.3). 

5.2.2 The isotropic background 

There are a number of other backgrounds to consider, such as cosmic rays or starburst 
galaxies. Since we expect all contributions of the EGB to add up to the experimentally 
observed flux, we must convolve the distributions above with the P{F ) of these other isotropic 
backgrounds. However, a complete model of these backgrounds is beyond the scope of this 
exploratory analysis. The one-point function of these other isotropic components is (for 
convenience) assumed Gaussian with a small, arbitrarily chosen, but non-negligible variance 
(m/o- r^j 10 3 ) and mean determined by requiring that the mean intensity due to -Pegb(T) (after 
adding the dark matter and blazar components) be equal to the experimentally determined 
value 6 x 10 -10 cm~ 2 s _1 sr -1 MeV -1 [65]. This Gaussian could be thought of as the central 
limit theorem approximation to the flux distribution of a large number of sources, since [i/a 
scales like y/N. 

We can also consider not convolving a dark matter component into the gamnra-ray 
background, to study (in Sec. 5.3) the distinguishability of an alternate hypothesis model 
with dark matter, P\it(F), from a null hypothesis model without dark matter P^ u n(F) 
(Both of these flux distributions are represented in Fig. 6). We then have 

TeGb(T) = PBlazar(T) *^ res t(E) * Pdu(F) , (5.6) 

which can be solved using the blazar P(F) model [Eq. (5.5)] above, as 

Pegb(F) = [e- 2 -°5(F) + (1 - e~ 2 -°)P™g ai (F)] *&est (F) * Pem(F) , (5.7) 

where the analytical factor of (1 — e^ 2 0 ) captures the normalisation of the Monte Carlo. It 
shall be crucial in what follows that the mean of Pegb(F) has been fixed as an experimental 
input, rather than allowed to vary freely as a model parameter. 

5.2.3 The combined astrophysical backgrounds 

The purpose of this section is to describe the influence of modelling choices on the shape of 
the flux distribution of the total unresolved EGB, as represented in Fig. 6. 

The thin peaks and power-law-tailed peaks at low and high flux in Fig. 6, correspond 
respectively to whether or not a blazar is present in the associated pixel. The relative peak 
heights are determined (to a first approximation) by the number of unresolved blazars per 
pixel, via the mixture coefficient e~ 20 . Broadening the peaks (by convolving the astrophysical 
components with a dark matter component) introduces a correction to this determination 
of height. The location of the mean of the low-flux peak is set via (i) Eqn. (5.6), (ii) the 
(assumed) mean of the blazar flux distribution (cf. Table 3), and (iii) the experimentally 
determined mean of the total flux distribution [65]. The difference between the dark matter 
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Figure 6: Predicted flux distribution Pegb(F) of the extragalactic ganrma-ray background, 
with (black) and without (green) a contribution from dark matter annihilations. The distri¬ 
butions have two peaks, based on whether or not a blazar is present in the associated pixel. 
The mean EGB derived from Fermi [65] is represented by the vertical line (red, dashed). A 
cross-section twice the canonical value was used to visually enhance the differences between 
these distributions. 


distribution’s mean and most probable fluxes (cf. Table 2) introduces an additional percent- 
level shift between the locations of the low-flux peaks with and without dark matter. 

The width and depth of the ‘gap’ between the high and low flux peaks is related, in 
the absence of dark matter, to the lower flux limit on blazars assumed in Sec. 5.2.1. This 
relation relies on the interplay between two effects: on one hand, extrapolating to lower 
fluxes increases the mean blazar contribution, shifting the location of the low-flux peak by an 
equal and opposite amount to even lower fluxes and widening the gap. On the other hand, 
extrapolating the faintest blazar contribution to lower fluxes also fills the gap with the fluxes 
from these faint blazars. The latter effect more than compensates for the former, such that 
overall the gap closes as it shifts to lower fluxes with an increasing blazar contribution. 

However, the contribution of the dark matter component to this gap is just as dramatic. 
Since the distribution is skewed, it broadens both peaks preferentially to higher fluxes; how¬ 
ever, the location of the low-flux peak is fixed by the constraint imposed on the mean of 
Pegb(F), while the high-flux peak is free to shift under this broadening. This widens the 
gap with increasing brightness of the dark matter annihilation signal. 

In addition to these shifts, the slope of the flux-tail due to dark matter annihilation 
(—2.5) is quite different from the one for blazars (—1.64); if the former component gives 
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a significant contribution to the gamma-ray background, this will deform the shape of the 
total distribution’s tail. The effect of changing the variance of the isotropic Gaussian, being 
of less interest than actually modelling these components, has not been investigated in this 
exploratory analysis. See Sec. 5.4. 

Note that all these features of the predicted gamma-ray flux distribution, both with 
and without dark matter, fit just on the edge of the experimental uncertainty on this mean 
value (6± 1 x 10 -10 cm -2 s -1 sr -1 MeV -1 , Ref. [65]), and may be open to exploitation by 
one-point function methods. 

As the sum of a diffuse background and a few unresolved sources, the gamma-ray back¬ 
ground in Eq. (5.6) has the same origin as the dark matter background. By adopting a 
wider variance for the Gaussian in Eq. (5.6), even the ‘Gaussian with a power-law tail’ form 
of the dark matter component can even be reproduced, corroborating our interpretation of 
these features of the dark matter P(F) and justifying our faint/bright E* analysis. One 
might worry that this ‘truly diffuse plus unresolved point sources’ form also undermines the 
prospects of an unambiguous dark matter detection above such a background. Indeed, model¬ 
fitting one-point techniques such as that described in Ref. [24] cannot separate degenerate 
models, especially given the angular resolution limitation (illustrated in Fig. 3) that prevents 
the dark matter E(E) power-law tail from contributing significantly to the observed flux. 
However, one need not worry too much: the asymmetry-induced shift of the peaks of E(E) 
may be a sufficiently distinctive feature to extract a dark matter signal from the Fermi data 
nonetheless. 

5.3 The photon count distribution P(C) 

The observable given by Fermi is not the gamma-ray flux E, but the discrete number of photon 
counts per pixel C. Photon arrival may then be modelled as a Poisson rate with a mean 
determined by the differential gamma-ray flux and the exposure e = (time) x (detector area) x 
(photon energy). For a five-year Fermi mission, correcting for the field of view, we have an 
exposure of e ~ 2.83 x 10 14 cm 2 s MeV sr pixel -1 . Marginalising over the uncertain flux 
distribution then gives 

P(C) = J P EGB (F)V(C\eF)dF. (5.8) 

This Poisson arrival uncertainty substantially smooths away the differences between the null 
and alternate flux models, as evidenced by Fig. 7. However, the percent-level shift between 
the low-flux peaks due to the dark matter distribution’s skewness survives, since a percent 
difference with C ~ 0(100) is still a few photons. There is also a larger, opposite shift in the 
point-source-driven high-flux tail due to our imposed value of the distribution’s mean. 

We can define, given our number of pixels IVpi x , the test statistic [25] 

2 _ (-^pix [EnuIi(C) — T’Alt(C')] 

X _ C V V^pi^NulKC) 

The choice of bounds for the sum over count bins is somewhat arbitrary; a formal optimisation 
of this test statistic would be beyond the scope of this analysis. We choose to focus on the 
peak 65 < C < 165 of the distribution, in which we anticipate sufficiently many pixels per 
‘count bin’ to trust a y 2 test. The lower panel of Fig. 7 illustrates the terms of the sum in 
this test statistic, from which one may obtain the p-value at which data with dark matter 
following Pak^C) exactly would reject our dark-matter-free null hypothesis. 
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Photon Count C [pixel 3 ] 

Figure 7 : Predicted count distributions of EBG photons with (black) and without (green) 
a dark matter component. The green bands represent the Poisson errors a oc y/N P(C) on 
the dark-matter-free model. The lower panel shows difference between the two models. A 
cross-section twice the canonical value was used to visually enhance the differences between 
these distributions. 


For the fiducial dark matter model, the percent-level shift between these peaks is just 
small enough that the null cannot be rejected by the data. Since the dark matter cross section 
(av) enters in our model only as a proportionality factor for the dark matter halo luminosity, 
we can rescale our Pdm(F) to inexpensively repeat this detectability study for higher values 
of the cross section. We can then forecast for which values of the cross section a dark matter 
component would become distinguishable from the background using the one-point function 
alone. This is summarised in Fig. 8, which shows that (given a perfect understanding of the 
backgrounds) the one-point function could probe a dark matter annihilation signal with a 
cross section roughly a factor two times larger than the canonical value 3 x 1CF 26 cm 3 s -1 . 
Fig. 8 also shows that our fiducial choice of lower blazar flux extrapolation is fortuitously 
close to a detectability optimum, but our forecast does not deteriorate much upon rescaling 
this value. 

This result, complementary to two-point function analyses, could even be strengthened 
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Cross Section {av) / 3 x 10 26 cm 3 s 1 

Figure 8: Predicted statistical significance of a (hypothetical) one-point-function-only de¬ 
tection of a dark matter annihilation signal above perfectly characterised astrophysical back¬ 
grounds, as a function of the dark matter cross section. Curves are labelled by the flux S m i n 
down to which the blazar distribution is extrapolated (see Table 3). Horizontal lines (blue, 
dashed) represent some common choices of confidence level. Including the energy-dependence 
of the flux distributions would improve these results, at the cost of a greater dependence on 
the annihilation spectrum. 


by including the energy dependence of the differential flux to break the degeneracy with the 
astrophysical backgrounds [28]. Such a study would remain sensitive to (but would allow a 
quantitative analysis of) the assumptions and uncertainties of the astrophysical background 
model. Yet, even without this spectral input, our forecasted one-point upper limit on the 
cross-section is on par with the most recent (spectral) constraints [5] based on the mean value 
alone. 

In addition to the extragalactic dark matter flux, there will be a component due to 
Galactic substructures. The one-point distribution of such a Galactic component has been 
predicted [25], and similarly features a power-law high-flux tail. Due to the energy spec¬ 
trum Eq. (3.8), if the mean intensity from subhalos at the anticenter integrated above 
10 GeV is ~ 10 -10 cnr 2 s~ 2 sr _1 [25], then the mean differential intensity at 1 GeV is 
(I) 10 12 cm 2 s 1 sr 1 MeV . This is of the same order of magnitude as the extragalac¬ 

tic component discussed above, and with the same high-flux F~ 2 ' 5 power-law tail. Thus, 
including the Galactic component would further enhance the expected signal-to-noise for po¬ 
tential detection. We finally note that the Galactic component will show a dipole feature, 
with more flux from the Galactic center than the anticenter, which can in principle be used 
to discriminate it from the isotropic extragalactic component. 

5.4 Caveats 

There are a number of caveats on the results presented in this study. Firstly, a large number 
of assumptions were used to simplify the hierarchical model without a proper sensitivity 
analysis. We assumed there is no scatter in the halo parameters [Eq. (2.6)], or any uncertainty 
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on the average number N' of halos in each pixel. An NFW profile was assumed for the 
halos despite the fact that Einasto and uncusped profiles would probably give less flux. The 
uncertainty incurred by extrapolating the mass function dN/dM and the boost models B(M ) 
down to M m i n , is compounded by our ignorance of the actual value of M m [ n . 

Secondly, we have not studied how our results depend on pixel size, particularly the 
effects of source extension. We have merely assumed (Sec. 2.2) that all sources are point-like, 
since there are on average only 0.28 extended dark matter sources per pixel, a negligible 
fraction of all N' = 7 x 10 21 halos (Sec. 3.1.1). However, extended sources must be either 
massive or nearby (Sec. 4.1), and therefore would tend to have large fluxes, affecting the 
distribution p£ >f *(F). Our point-source-based Pi(F) is therefore not applicable to these 
objects at high flux. That this compromises the analysis should be obvious from Fig. 4: the 
one-point functions do not account for the clusters’ extension, and the clusters’ fluxes clearly 
do not live in the domain of the PDF they should be drawn from. Thus, extended sources 
should be dealt with in a complete analysis. Some of the elements of such an analysis (such 
as an energy-dependent angular resolution, a substructure-boost-dependent mass-to-solid- 
angle conversion, or a redshift-dependent mass threshold M^xtiz)) have been presented in 
the main text. Note, however, that the distribution of faint and distant sources P f<f *(F ) 
used to derive Qq^y*{E) is not compromised by the extension of bright and nearby sources; 
as long as the faint/bright split in the method of Sec. 4.2.1 excludes faint extended sources 
(which can be guaranteed by a suitable choice of F 1 *), only the Monte Carlo based on bright 
sources would need be revised. 

Thirdly, we have not considered the energy spectrum of our dark matter annihilation 
signal, besides noting that it is a very relevant quantity in Fig. 1 (which is clearly insufficient). 
Using differential fluxes throughout this study is a first step in this direction. On one hand, 
the choice of the particle dark matter spectrum dN/dE is mostly unrelated to structure 
formation and only contributes an energy-dependent normalisation which, from a particle 
physics perspective, is almost completely arbitrary. Hence, it would not have much predictive 
power and matter in practice only for data fitting. On the other hand, at higher energies 
the gamma-ray absorption matters, and a more sophisticated model than Eq. (3.9) for this 
quantity would be required. The energy resolution of the instrument (A E/E = 9% at 
1 GeV [1]) would also need to be accounted for in the model. 

Finally, our Gaussian model for the non-blazar isotropic components of the EGB is 
clearly inadequate (and even the blazar model is somewhat simplistic), since rigorously ac¬ 
counting for all the astrophysics would require an entirely separate analysis. Consequently 
the one-point functions and forecasted limits that depend on this input must be understood 
as exploratory and methodologically illustrative. 

6 Conclusions 

We constructed a hierarchical model that predicts, using analytical models of ACDM struc¬ 
ture formation, the flux distribution of gamma rays from extragalactic dark matter annihi¬ 
lation in unresolved point sources. The uncertainties on this flux subject to the modeling 
choices we studied are typically percent-level; in the case of the substructure boost function, 
they remain smaller than a factor of three. We then compute, without requiring any addi¬ 
tional physical assumptions, the flux distribution per pixel P{F), which has the characteristic 
form of an isotropic diffuse Gaussian matched at high flux to the point-source distribution 
with a power-law slope of —2.5. This distribution is non-Gaussian and asymmetric; however 
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the most likely flux and the mean flux are comparable at the percent-level in all but the 
optimistic boost model, salvaging previous ‘mean intensity’ constraints on the dark matter 
properties from this potential systematic effect. 

The fluxes predicted for our fiducial model lie just within the reach of the Fermi-LAT, 
and should be observable by the tenth year of the mission. We also showed that the dis¬ 
tinctive features of the power-law-tailed Gaussian distribution all live above Fermi’s angular 
resolution. Therefore, the extragalactic gamma-ray emission due to dark matter annihilation 
constitutes an irreducible and significant background for point-source annihilation searches 
with clusters or dwarf spheroidals. Ironically, an optimistic boost model would be detrimen¬ 
tal to these searches, by deteriorating the signal-to-noise of these point sources (to unity or 
worse for galaxy clusters). 

We also discussed the astrophysical backgrounds from which a dark matter annihila¬ 
tion signal would need to be extracted. These include unresolved blazars (which contribute 
an order of magnitude more flux than the fiducial dark matter model) and other diffuse 
components, which were all convolved together into a total model for the gamma-ray back¬ 
ground. The scarcity of unresolved blazars make this distribution quite rich in features; most 
prominently, it has two distinct peaks of most probable fluxes, the inter-peak gap being very 
sensitive to the dark matter component. 

Even accounting for the Poisson noise of photon arrivals that come with such low fluxes, 
a contribution to the gamma-ray background of the order of a vanilla WIMP model may be de¬ 
tectable above well-characterised astrophysical backgrounds using the flux distribution alone. 
Using the energy-dependence of the flux distribution should further break the degeneracy be¬ 
tween the components of the gamma-ray background, and should allow one-point function 
methods to complement and strengthen existing constraints set by two-point-function anal¬ 
yses. 
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A Failure of the central limit theorem 

When constructing P(F) from many individual sources (each an independent and identically 
distributed realisation of P\(F)), it is easy to show that we expect significant deviations from 
a Normal distribution, even for the large number N' ~ 7 x 10 21 of sources per pixel. The 
ratio between the mean flux and the standard deviation of the fiducial model is approximately 
( ii/ct)n’ = 0.23, which (using the Gaussian’s associated cumulative distribution) would firmly 
place 40% of a Gaussian P(F) at negative fluxes. This is not only mathematically impossible 
(given that the sum of positive random variables must be positive), it is also physically 
nonsensical. Note that using distributions with widths instead of delta functions for P(N') 
or V ( k | N’ ) makes this problem even worse by broadening the unphysical ‘Gaussian’ Pf.. 

We can also show that we expect deviations from a Gaussian with slightly more rigour: 
The /c-autoconvolution definition of Pk{F) = (Pi)* k (F) prompts us to work in the Fourier 
space of probability distributions, i.e., with the characteristic functions and 4>k = (4>i) k ■ 
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Marginalising Pp.(P) with the Poisson distribution V (k\N') gives [30]: 


P(F\N') = £ 


00 „-N' jy/k 
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k =0 
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where we have used linearity of the (inverse) Fourier transform and have recognised the power 
series expansion of the exponential function. By Taylor expanding the characteristic function 
4>i, we generate (by construction) the first moments of Pi: 

= 1 +E(Pi)(it) + V(Pi)^ + o(-it 3 ), (A.2) 

such that the characteristic function associated to P(F |/r) becomes (to second order in t) 

<t>p(F\ N >) = (Ve(Pi)] (it) - [JV'V(Pi)] e -/v'o(-it 3 ) (A.3) 

We recognise this first exponential as the characteristic function of a Normal distribution, so 


P(F\N') » Q (F\N'E(Pi), N'V(Pi)) *P _1 


e N'o(-it 3 ) 


(A.4) 


If P(N') is a (5-function we recover upon marginalisation the near-Gaussian form of P(F) 
anticipated in the main text, with explicit deviations due to higher moments convolved in. 
Proving that these deviations lead to a power-law tail at high flux is of less interest than 
the observation, in the main text, that in this regime the flux is almost completely due to a 
single bright source. This leads, in the next section, to a Monte Monte CarloCarlo analysis 
of these corrections. 

The analysis above does not assume N',k —> 00 . This reflects the fact that we have not 
derived the Central limit Theorem (CLT), the ‘Gaussian * Corrections’ form was derived for 
a finite number of sources. However, the CLT is still valid asymptotically as k —> 00 since 
Pi does have a finite variance despite its power-law-like behaviour. 


B Our method 

B.l Rationale of our method 

The deviations from Gaussianity in the sum P*. of random variables are bounded by the 
Berry-Esseen (BE) theorem [31, 63], which relates the Kolmogorov-Smirnov distance (the 
largest deviation at any point), to the ‘absolute skewness’ 5 of P\ (F) and the large but finite 
number of dark matter halos k: 

sup |CDFfc(F) — < h(P)| < —v=- , Pl = E{lX ~/ l[ '\ (B.l) 

F Vk cr 1 

5 Formally, pi is the sum of the partial skewnesses above and below the mean. We simply call pi the 
skewness. 
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where C ~ 0.5, and where CDF^. and &(F) are the cumulative distribution for Pf. and the 
Normal distribution, respectively. We note that the absolute value in p\ distinguishes it from 
the conventional skewness 71 , and gives the inequality p\ > 71 . For our power-law-like Pi ( F ), 
the skewness is huge and the BE bound Cp\k~ l t 2 is uninformative: the Kolmogorov-Smirnov 
distance (a quantity definitionally less than one), is ‘constrained’ by the first few moments 
of our fiducial P\ (F) to be less that about 10 8 . 

This perspective suggests splitting P\(F) into low-flux and high-flux contributions, to 
reduce the skewness of each contribution. In the following, we (i) derive and (ii) cross-check 
the Monte Carlo method presented in the main text. 


B.2 Deriving the form of Pk(F ) 

We derive the behaviour for P as follows: marginalise Pi (F) into high and low flux contri¬ 
butions, 

Pi(F) = (1 - e)Pf <F *(F) + eP F>F * (F ), (B. 2 ) 

where P 1 *’ *(F) are the normalised distributions of fluxes below and above the trun¬ 

cation flux jF* (see Fig. 2 of the main text), and the fraction of high-flux sources, 

ro o 

e = e(F*) = / Pi{F)dF, (B.3) 

JF, 

is the relative normalisation of these two distributions. Clearly the smaller we choose I 7 *, the 
larger e becomes. The number of faint sources A;* is then simply A:(l — e). 

Pk is the /c-autoconvolution of the sum above: we may then write Pk as the sum 


Pk = 


k\ 


k 


k—i 


pF<F * p_F>F* 

F-i * F 


(B.4) 


This can easily be shown by induction or by using the binomial theorem on the associated 
characteristic function 


fa = ( fa) k = ((1 - + ecj) F>F *) k . 


(B.5) 


Unless k S> i, the contribution of the zth term is vanishingly small. Physically, this 
reflects the low probability of having many sources from the high-flux tail of the ‘power-law’ 
Pi(F). The low-flux autoconvolution P^ F * is very nearly Gaussian when F M is small (i.e. 
when the variance of P F<F * is small). Each term of the sum is then an isotropic background 
Gaussian, convolved with high fluxes from point sources (drawn from P F>F *). But, since 
k 3> i the low-flux contribution P F -^ F * factorizes (see Eqn. (B.7) below). 

How many terms do we need to consider for the gargantuan sum [Eq. (B.4)] to converge? 
We expect a number ke of high flux sources on average. In other words, the number i of high 
flux sources in a pixel is drawn from the Poisson distribution V (i\ke). This is reflected by 
the binomial expansion: for 1 e (which also implies k 3 > i) we have 


A ; 1 . 

— (l-e) fc -* e l 


(key 


= e kt V (i\ke ), 


i!(k-i)r~ i\ 

such that (reassembling all the elements of the discussion above) 
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p ~ pF<F t ke 

Fk ~ F, ★ e 
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(B. 6 ) 


(B.7) 
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Therefore, the physically relevant terms in the sum are those in the 5a band ke — 5 V~ke < 
i < ke + 5 Vke of the distribution V (i\ke). 

The bracketed sum, even with its large number k of terms, can be evaluated by Monte 
Carlo: for each realisation, draw i from a Poisson distribution, then draw that many fluxes 
from P k >F *, and sum over the fluxes. The histogram of F^ F * over many realisations then 
approximates the bracketed term (i is marginalised automatically by the algorithm). We 
then rescale by a constant e ke (absorbed into the normalisation of Pk(F)) and convolve with 
the isotropic background of low-flux sources (just a thin delta-function in practice) to get P^. 

B.3 Cross checks 

We should, of course, check that our choice of F* is physically sensible. If we want to split 
Pk{F) into an isotropic background and candidate point sources, we must verify that our 
faint sources do not contribute more than zero or one photons each. Recall that a single 
GeV photon per pixel in 5 years of Fermi data corresponds to a differential flux/intensity of 
F = 6 X 10 -6 m -2 yr -1 MeV -1 pixel -1 . Choosing F* « l.lx 10 -8 m -2 yr -1 MeV -1 pixel -1 , 
over two orders of magnitude smaller, guarantees that faint sources are indeed faint. This 
value of F 1 * was in fact chosen algorithmically, such that the BE bound for P f<f *(F) be 
Cpk~ 1 / 2 = 0.005. The applicability of the CLT in this low-flux component is thereby justified. 

We must also check that our choice of F 1 * does not noticeably influence the final distribu¬ 
tion P{F ): This cut is a mathematical artifact. Nature does not fundamentally distinguish 
halos that appear ‘faint’ or ‘bright’ to a specific observer, with an arbitrary F 1 *. In Fig. 9 
we present a (local) sensitivity analysis using three different F 1 * cuts, which shows that even 
a cutoff two orders of magnitude greater than ours (and just below the one-photon-per-five- 
years threshold) would have given satisfactory results (although the BE bound on the low-flux 
dark matter component would have been appreciably weaker). The relative error near the 
peak of the distribution is percent-level despite order-of-magnitude variations, even for a rel¬ 
atively small Monte Carlo (10 6 realisations). This error has an expectation consistent with 
zero and scales with the size of the Monte Carlo, as expected of an Monte Carlo sampling 
noise. 

Although P(F) is increasingly poorly-sampled by our Monte Carlo (Fig. 9) as the flux 
increases, the behaviour of this tail remains well-known: since the probably of a high-flux 
source is relatively low, the tail of P(F) should look like P±(F). This behaviour is illustrated 
in Fig. 10. The power-law tail continues until its cutoff (Fig. 1), so the power-law tail of 
P(F) spans roughly ten orders of magnitude in flux. In light of the Gaussian intuition of 
exponential tails, this power-law tail explains how our P(F) can support such a large variance 
despite the sharpness of its peak. 

C Angular resolution limit 

We can compute the probability of seeing a dark matter signal above a given flux, by looking 
at the exceedance (complementary cumulative) distribution ’F(F’) associated to P(F): 

roo 

^(Fjnax) = / P{F)dF. 

•' fmax 

The probability of not realising this high flux tail in Np- lxe i trials is then given by the binomial 
distribution, and we want to solve the following for F max : 

I3(0|IV pixel ,T(F max )) < a , a = {0.05,0.01,0.001}. (C.l) 
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Figure 9: Local sensitivity analysis of our Monte Carlo on the cutoff flux F*, with 10 6 
samples of P F>F * for each F*. The values of the cutoff F* quoted in the legend are in units 
of nr -2 yr -1 MeV -1 pixel -1 . The error in the lower panel is relative to the average of the 
three Monte Carlo simulations. 


For k = 0 successes and p = 'L(F max ) < 1, we can expand this in a binomial series that we 
can truncate at first order when T(F max ) IVF 1 ^: 

1 - AWci^max) + O(fVpi X elT(F max )) 2 < a. (C.2) 

When F max lies in the power-law tail of F(F), we can use the high-flux-tail equivalence 
(Fig. 10) of P(F) and Fl(F) to find, approximately, 

1 - a 

/ Pi(F)dF< — . (C.3) 

JF max JV pixel 

For a power-law-like Fi(F) ~ AF^,^ < 0 this integration yields 

roo 

/ AF^dF = -^(F max )^+V(7 + 1) = -PmaxFi(F max )/( 7 + 1). (C.4) 

J ^max 
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Intensity F/r2pi x [cm 2 s 1 sr 1 MeV x ] 

Figure 10: Unnormalised P(F ) and P\ (F) (solid, black) for the fiducial model. The scaling 
factor to achieve the matching between the two distributions is not a free parameter, it is 
the number of halos k — fc* = 1864 contributing to the high-flux tail in the fiducial model. 
We show also the (unnormalised) Gaussian with equivalent mean and variance predicted by 
a naive CLT (green, dashed). 


In the limit that (a, 7 ) —> (0, —2.5) we reproduce Eq. (4.4) from the main text. More 
generally, for a —> 0 and 7 = —2.5 + dy, we have 

Tmax-Pi (F max ) < - -p- . (C.5) 

1 ''pixel pixel 

Since (as visible in Fig 1) the highest fluxes in the tail of Pi(F) have a steeper log-slope than 
—2.5 (i.e. 5y < 0), the angular resolution limit of Eq. (4.4) is actually more optimistic than 
would be warranted by a more precise calculation. 
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